# -*- coding: utf-8 -*-
"""
Created on Wed Jun 05 14:44:21 2013

@author: Maxim
"""

import propeller
import engine
from scipy.optimize import fsolve,minimize_scalar,minimize
from numpy import zeros,linspace,array,array_equal,transpose,vstack
import matplotlib.pyplot as plt


class Thrust:
    def __init__(self,propeller,engine):
        self.prop = propeller
        self.prop.neval = 0
        self.engine = engine
        self.S = 11.4
        self.W = 620.0*9.81
        self.CD0 = 0.04
        self.k = 0.01
        self.beta = propeller.beta75
        self.nAnalysis = 0
    def set_pitch_angle(self,beta):
        self.prop.set_beta(beta)
        self.beta = beta
    def _vel_obj(self,x,pset,rho):
        rpm, V = x[0],x[1]
        beta = self.beta
        data = self.prop.analyze_prop(beta,rpm,V,rho)
        Tav = data[0]
        Preq = data[1]
        Pav = self.engine.get_availablePower(rho,rpm,pset)
        q = rho*V*V/2.0
        qS = q*self.S
        CL = self.W /qS
        CD = self.CD0 + self.k*CL*CL
        Treq = CD*qS
        return [Tav-Treq,Pav-Preq]
    def velocity_at_power_setting(self,pset,rho):
        x = fsolve(self._vel_obj,[1000.0,40.0],args=(pset,rho),xtol=1e-8)
        return x
    def max_efficiency(self,pset,rho):
        def obj(beta,pset,rho):
            self.nAnalysis += 1
            self.set_pitch_angle(beta)
            x = self.velocity_at_power_setting(pset,rho)
            rpm, V = x[0],x[1]
            data = self.prop.analyze_prop(beta,rpm,V,rho)
            return -data[-1]
        betaBest = minimize_scalar(obj,bounds=(15.0,40.0),args=(pset,rho),method='bounded')
        self.set_pitch_angle(betaBest.x)
        rslt = self.velocity_at_power_setting(pset,rho)
        print rslt
        print betaBest
        return rslt
    
    def max_efficiency2(self,pset,rho):
        def obj(x,pset,rho):
            self.nAnalysis += 1
            beta,rpm,V = x[0],x[1],x[2]
            data = self.prop.analyze_prop(beta,rpm,V,rho)
            return -data[-1]
        
        def constr(x,pset,rho):
            x2 = [x[1],x[2]]
            self.set_pitch_angle(x[0])
            d = self._vel_obj(x2,pset,rho)
            return d
        x0 = [20.0,1000.0,30.0]
        bnds = ((15.,40.),(500.,2300.),(25.,65.))
        cnstr = {'type':'eq','fun':constr,'args':(pset,rho)}
        x = minimize(obj,x0,method='SLSQP',bounds=bnds,constraints=cnstr,
                     args=(pset,rho))
        print x
        print constr(x.x,pset,rho)


class FixedPitchThrust:
    def __init__(self,engine):
        self.engine = engine
        self.prop = None
        self.thicknessOld = zeros(3)
        self.S = 11.4
        self.W = 620.0*9.81
        self.CD0 = 0.04
        self.k = 0.03985
    def _vel_obj(self,x,pset,rho,weight):
        rpm, V = x[0],x[1]
        beta = self.beta
        data = self.prop.analyze_prop(beta,rpm,V,rho)
        Tav = data[0]
        Preq = data[1]
        Pav = self.engine.get_availablePower(rho,rpm,pset)
        q = rho*V*V/2.0
        qS = q*self.S
        CL = weight /qS
        CD = self.CD0 + self.k*CL*CL
        Treq = CD*qS
        return [Tav-Treq,Pav-Preq]
    def _vel_obj2(self,rpm,V,pset,rho,weight):
        beta = self.beta
        data = self.prop.analyze_prop(beta,rpm,V,rho)
        Tav = data[0]
        Preq = data[1]
        Pav = self.engine.get_availablePower(rho,rpm,pset)
        q = rho*V*V/2.0
        qS = q*self.S
        CL = weight /qS
        CD = self.CD0 + self.k*CL*CL
        Treq = CD*qS
        return ((Tav-Treq)**2.0+(Pav-Preq)**2.0)**.5
    def thrust_at_velocity(self,V,pset,rho,weight):
        rpm = fsolve(self._vel_obj2,[1000.0],args=(V,pset,rho,weight),xtol=1e-10)
        return rpm
    def velocity_at_power_setting(self,pset,rho,weight):
        x = fsolve(self._vel_obj,[1000.0,40.0],args=(pset,rho,weight),xtol=1e-10)
        rpm, V = x[0],x[1]
        return rpm,V


def run_test1():
    pset = 50.0
    prop = propeller.load('neuform')
    eng = engine.EngineExperimental('ROTAX_912')
    thrust = Thrust(prop,eng)
    thrust.set_pitch_angle(20.0)
    x = thrust.velocity_at_power_setting(pset,1.05)
    lb = eng.get_rpm_min(pset)
    ub = eng.get_rpm_max(pset)
    rpm = linspace(lb,ub,10)
    power = list()
    for r in rpm:
        power.append(eng.get_availablePower(r,pset,1.05))
    plt.figure(1)
    plt.grid(True)
    plt.hold(True)
    plt.plot(rpm,power,'bo')
    plt.plot(x[0],eng.get_availablePower(x[0],pset,1.05),'ro')
    plt.show()

def run_test2():
    eng = engine.load('ROTAX_912ULS')
    prop = propeller.load('neuform')
    thrust = Thrust(prop,eng)
    pset = 95.0
    rho = 1.05
    prop.set_beta(15.0)
    print 'direct'
    thrust.max_efficiency(pset,rho)
    print thrust.prop.neval
    thrust.prop.neval = 0
    print '\n IDF'
    thrust.max_efficiency2(pset,rho)
    print thrust.prop.neval
    print eng.get_fuelFlow(rho,2300.,pset)


def test_03():
    eng = engine.load('ROTAX_912ULS')
    prop = propeller.load('MDO_optimum')
    pset = 100.0
    rho = 1.2255
    weight = 620.0*9.81
    thrust = FixedPitchThrust(eng)
    thrust.prop = prop
    thrust.beta = 28.5
    rpm = thrust.thrust_at_velocity(50.0,pset,rho,weight)
    print rpm
    
if __name__=="__main__":
    run_test2()

#        return T,P,Jreq/Jcorr,CT,CP,effy
